
%%% The Ambiguity of Fishing for Fun, October 2022
%%% This function evaluates the loglikelihood and its gradient for the alpha-maxmin utility model

function[LnL , gradient]=likelihoodGrad_ambiguity(R,Z,Y,C,I,IJ,ID,parameters)

%Calculating min and max number of kept and released fish for summer flounder
max_G=R(:,3);
min_G=R(:,4);
max_Q=R(:,1)-R(:,3);
min_Q=R(:,2)-R(:,4);

delta=[parameters(3), parameters(4)]'; 
vartheta=[parameters(8),parameters(9),parameters(10),parameters(11),parameters(12),parameters(13), parameters(14)];

W_min=parameters(1).*min_Q+parameters(2).*min_G+Z*delta;
W_max=parameters(1).*max_Q+parameters(2).*max_G+Z*delta;
u_min=(1-exp(-parameters(5).*W_min))/parameters(5);
u_max=(1-exp(-parameters(5).*W_max))/parameters(5);

F=parameters(15).*u_min+(1-parameters(15)).*u_max+parameters(6).*C; 
EF=exp(F);

%Partial derivatives with respect to theta, r (risk aversion) and alpha (ambiguity aversion)
f_theta1=parameters(15).*(1-parameters(5)*u_min).*min_Q+(1-parameters(15)).*(1-parameters(5)*u_max).*max_Q;
f_theta2=parameters(15).*(1-parameters(5)*u_min).*min_G+(1-parameters(15)).*(1-parameters(5)*u_max).*max_G;
f_theta3=Z(:,1).*(parameters(15).*(1-parameters(5)*u_min)+(1-parameters(15)).*(1-parameters(5)*u_max));
f_theta4=Z(:,2).*(parameters(15).*(1-parameters(5)*u_min)+(1-parameters(15)).*(1-parameters(5)*u_max));
f_theta5=parameters(15).*((((1/parameters(5))*(W_min+1/parameters(5))).*(1-parameters(5)*u_min))-1/(parameters(5)^2))+(1-parameters(15)).*((((1/parameters(5))*(W_max+1/parameters(5))).*(1-parameters(5)*u_max))-1/(parameters(5)^2));
f_theta15=u_min-u_max; 

%Calculating the opt-out utility
Eopt_out=exp(parameters(7).*ones(length(C),1)+Y*vartheta'+parameters(6).*C);

%Separating trips A and B from the opt-out trip alternative
ef=EF.*(1-IJ);
eopt_out=Eopt_out.*IJ;
Denom=ef+eopt_out;
InactiveCosts=C.*IJ;
sumtheta1=f_theta1.*EF.*(1-IJ);
sumtheta2=f_theta2.*EF.*(1-IJ);
sumtheta3=f_theta3.*EF.*(1-IJ);
sumtheta4=f_theta4.*EF.*(1-IJ);
sumtheta5=f_theta5.*EF.*(1-IJ);
sumtheta6=C.*EF.*(1-IJ);
sumtheta15=f_theta15.*EF.*(1-IJ);

%Forming the denominator
sume=accumarray(ID, Denom);
sumeInactiveCost=accumarray(ID, InactiveCosts);
sumActiveTrips=accumarray(ID, ef);
sumInactiveTrips=accumarray(ID,eopt_out);
Sumtheta1=accumarray(ID,sumtheta1);
Sumtheta2=accumarray(ID,sumtheta2);
Sumtheta3=accumarray(ID,sumtheta3);
Sumtheta4=accumarray(ID,sumtheta4);
Sumtheta5=accumarray(ID,sumtheta5);
Sumtheta6=accumarray(ID,sumtheta6);
Sumtheta15=accumarray(ID,sumtheta15);
SUME=kron(sume,[1 1 1]'); 
SUMEInactiveCost=kron(sumeInactiveCost,[1 1 1]');
SUMEActiveTrips=kron(sumActiveTrips,[1 1 1]');
SUMEInactiveTrips=kron(sumInactiveTrips,[1 1 1]');
SUMEtheta1=kron(Sumtheta1,[1 1 1]');
SUMEtheta2=kron(Sumtheta2,[1 1 1]');
SUMEtheta3=kron(Sumtheta3,[1 1 1]');
SUMEtheta4=kron(Sumtheta4,[1 1 1]');
SUMEtheta5=kron(Sumtheta5,[1 1 1]');
SUMEtheta6=kron(Sumtheta6,[1 1 1]');
SUMEtheta15=kron(Sumtheta15,[1 1 1]');

F=ef./SUME; 
FTerm=F(F~=0);
FTermI=I(F~=0);
S=eopt_out./SUME; 
STerm=S(S~=0);
STermI=I(S~=0);

%Calculating lnL 
LnL=-(sum(log(FTerm).*FTermI)+sum(log(STerm).*STermI));

%Calculating the gradient of lnL
if nargout>1
    gradient(1)=-(sum(I.*(1-IJ).*(f_theta1-(SUMEtheta1./SUME)))-sum(I.*IJ.*(SUMEtheta1./SUME)));
    gradient(2)=-(sum(I.*(1-IJ).*(f_theta2-(SUMEtheta2./SUME)))-sum(I.*IJ.*(SUMEtheta2./SUME))); 
    gradient(3)=-(sum(I.*(1-IJ).*(f_theta3-(SUMEtheta3./SUME)))-sum(I.*IJ.*(SUMEtheta3./SUME)));
    gradient(4)=-(sum(I.*(1-IJ).*(f_theta4-(SUMEtheta4./SUME)))-sum(I.*IJ.*(SUMEtheta4./SUME))); 
    gradient(5)=-(sum(I.*(1-IJ).*(f_theta5-(SUMEtheta5./SUME)))-sum(I.*IJ.*(SUMEtheta5./SUME)));
    gradient(6)=-(sum(I.*(1-IJ).*(C-(SUMEtheta6+SUMEInactiveCost.*SUMEInactiveTrips)./SUME))+sum(I.*IJ.*(SUMEInactiveCost-(SUMEtheta6+SUMEInactiveCost.*SUMEInactiveTrips)./SUME))); %Corresponding to parameters(6)
    gradient(7)=-(sum(I.*IJ.*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*(SUMEInactiveTrips./SUME)));
    gradient(8)=-(sum(I.*IJ.*Y(:,1).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,1).*(SUMEInactiveTrips./SUME))); 
    gradient(9)=-(sum(I.*IJ.*Y(:,2).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,2).*(SUMEInactiveTrips./SUME))); 
    gradient(10)=-(sum(I.*IJ.*Y(:,3).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,3).*(SUMEInactiveTrips./SUME))); 
    gradient(11)=-(sum(I.*IJ.*Y(:,4).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,4).*(SUMEInactiveTrips./SUME))); 
    gradient(12)=-(sum(I.*IJ.*Y(:,5).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,5).*(SUMEInactiveTrips./SUME))); 
    gradient(13)=-(sum(I.*IJ.*Y(:,6).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,6).*(SUMEInactiveTrips./SUME))); 
    gradient(14)=-(sum(I.*IJ.*Y(:,7).*(SUMEActiveTrips./SUME))-sum(I.*(1-IJ).*Y(:,7).*(SUMEInactiveTrips./SUME))); 
    gradient(15)=-(sum(I.*(1-IJ).*(f_theta15-(SUMEtheta15./SUME)))-sum(I.*IJ.*(SUMEtheta15./SUME))); 
end
